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f™^ Abstract. We derive energy- norm a posteriori error bounds for an Euler 

^_^ time-stepping method combined with various spatial discontinuous Galerkin 

f*^ schemes for linear parabolic problems. For accessibility, we address first the 

^sj spatially semidiscrete case, and then move to the fully discrete scheme by 

introducing the implicit Euler time-stepping. All results are presented in an 

13 abstract setting and then illustrated with particular applications. This enables 

^^H the error bounds to hold for a variety of discontinuous Galerkin methods, pro- 

| .^] vided that energy-norm a posteriori error bounds for the corresponding elliptic 

problem are available. To illustrate the method, we apply it to the interior 
penalty discontinuous Galerkin method, which requires the derivation of novel 
a posteriori error bounds. For the analysis of the time-dependent problems we 
^.^ use the elliptic reconstruction technique and we deal with the nonconforming 

^0^ part of the error by deriving appropriate computable a posteriori bounds for 

/'^ it. We illustrate the theory with a series of numerical experiments indicating 

jT , the reliability and efficiency of the derived a posteriori estimates. 

r^ Finite element, discontinuous Galerkin, error analysis, a posteriori, time depen- 
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dent problems, parabolic PDE's, upper bounds, nonconforming methods, time step- 
ping, Euler scheme 

~~ 1. Introduction 

%!^J Adaptive methods for partial differential equations (PDE's) of evolution type 

p,^ have become a staple in improving the efficiency in large scale computations. Since 

\^ the 1980's many adaptive methods have been increasingly based on a 'posteriori 

CNj error estimates, which provide a sound mathematical case for adaptive mesh refine- 

\l ment, which can be decomposed in spatially and temporally local error indicators. 

^^ In the context of parabolic equations, a posteriori error estimates have been derived 

^^ for various norms since early 1990's [TTIISS]- Inspired by the milestones set recently 

for the mathematical theory of convergence for adaptive methods in elliptic prob- 
lems [33 m [T3] , there has been a recent push for similar results for parabolic prob- 
lems calling to a closer understanding of a posteriori error estimates [15j |4TJ [71 [29l 
e.g.]. Most results in this area cover simple time-stepping schemes and a conforming 
space discretization. The extant literature on a posteriori error control for noncon- 
d forming spatial methods can be grouped in a handful of works [101 dSl HSl [Ml HI] • 

In [40] a posteriori L2(H )-norm error bounds for a spatially semidiscrete method 
via interior penalty discontinuous Galerkin (IPDG) methods are derived and used 
heuristically in the implementation of the fully discrete scheme. In [THlllS] L2(L2)- 
norm error bounds for IPDG are obtained using duality techniques, while in [51] . 
a posteriori error bounds are presented for a fully discrete method consisting of a 
backward Euler time-stepping and linear Crouzeix-Raviart elements in space. Note 
that none of the papers in the literature, to our knowledge, cover the case of a 
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posteriori energy-norm error bounds for fully discrete schemes with discontinuous 
Galerkin methods, which is the chief objective of our paper. 

Discontinuous Galerkin (DG) methods are an important family of nonconform- 
ing finite element methods for elliptic, parabolic and hyperbolic problems dating 
back to 1970's and early 1980's [35J |371 13 ilE]. DG methods have undergone 
substantial development in the recent years [TBJ SI [SHI EH [231 e-g-i and references 
therein]. The practical interest in DG methods owes to their flexibility in mesh 
design and adaptivity, in that they cover meshes with hanging nodes and/or locally 
varying polynomial degrees. DG methods are thus ideally suited for /ip-adaptivity 
and provide good local conservation properties of the state variable. Moreover, 
in DG methods the local elemental bases can be chosen freely for the absence 
of interelement continuity requirements, yielding very sparse — in many cases even 
diagonal — mass matrices even with high precision quadrature. Note also that DG 
methods are popular due to their very good stability properties in transport- or 
convection-dominated problems [16' ; the a posteriori error analysis of convection- 
dominated problems is, however, beyond the scope of our study and we concentrate 
on diffusion-only parabolic equations. 

Our main results are a posteriori error bounds in the energy norm for a family 
of fully discrete approximations of the following PDE problem — in ^we gather the 
functional analysis notation and background. 

1.1. Problem (linear parabolic boundary-initial value problem). Given an open 
(possibly curvilinear) polygonal domain Q C R'', d = 2, 3, a real number T > 0, 
two (generalized) functions 

/eLoo(0,r;L2(/?)) and a eL^i^ x {0,T)f'"', (1.1) 

such that a{x, t) is symmetric positive definite for almost all (x, i) G i? x [0, T], find 
a function u G L2(0,r;Hj(l2)), 

5tueLo,(0,T;H-i(l2)) (1.2) 

and such, that 

dtu-W ■ (aVu) = f on n X (0, T] , 
w(0) = Mo on n, and u\g^^ (t) = 0, for t e (0, T] . 

In [|2] we propose a class of numerical methods for solving this problem. These 
methods consist in a backward Euler time-stepping scheme in combination with 
various choices of spatial DG methods. Our emphasis is on the widely applied 
IPDG method [1IM1I21- 

We consider the notation of Problem [lT] to be valid throughout the paper and u 
denotes the solution of problem (1.3). Although the assumption / e Loo(0,T; L2(J7)) 
may be weakened — provided a posteriori error estimates for the corresponding spa- 
tial finite element method can be obtained for such weak data — we refrain from 
doing it for simplicity's sake. In fact, we consider / to be piecewise continuous in 
time with a finite number of time-discontinuities and with the implied constraints 



(1.3) 



on the time partition, to be discussed in ^2.4 The matrix-valued function a(-,i), 
for each t e (0, T) is allowed to have jump discontinuities; the set of spatial dis- 
continuities of a will be considered to be aligned with the finite element meshes. 
For simplicity, we shall assume that a is continuous in time, but a finite number 
of discontinuities can be accounted for easily, as long as these occur at the points 



of the time partition in the fully discrete scheme. The PDE (1.3) is assumed to be 



uniformly elliptic in the sense that the supremum and the infinium of the set 

[(a(x,t)C)-C/\Cf- CeR',{x,t)enx{Q,T]} (1.4) 
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are both positive real numbers. As for the boundary values, we remark that our 
approach can be appropriately modified in order to extend homogeneous to general 
time-dependent Dirichlet boun dary values. Under the assumptions made so far, 



we have that the solution to (1.3) exists and satisfies u e C(0,T; Hq(/?)) and 
5tueL2(0,T;L2(l2)) M- 

In line with a unified approach to a posteriori error analysis for elliptic-problem 
DG methods [U El US] our discussion will be presented first in an abstract setting. 
Our results are then shown to be applicable to a wide class of DG methods pro- 
vided that a posteriori error bounds for the corresponding steady-state problem are 
available. 

We stress that, although we focus on a posteriori error bounds for spatial DG 
methods, our abstract results can be applied to a wider class of nonconforming 
methods (other than DG methods), provided they satisfy certain requirements. 
More specifically, given a particular nonconforming finite element space Sh, assume 
that: 

(a) for each Z e Su it is possible to decompose it as 

Z = Zc + Zd such that Zc e Hj(i7) n Sh, (1.5) 

where Zc and Zd are called Z's conforming part and the nonconforming part, respec- 
tively. This decomposition is an analytic device and is not needed for computational 
purposes. The only requirement on this decomposition is the ability to quantify cer- 
tain norms of Z^ in terms of Z, as found in the literature [51 [551 [231 e.g.], as well 
as our Lemma 14.31 

(b) Given a function z € HJ(^), and let Z e S'h be the corresponding Ritz- 
projection via the finite element method, it is possible to bound the norm of the 
error Z — z, using a posteriori error estimators for the steady state problem. 

A key tool in our a posteriori error analysis is the elliptic reconstruction tech- 
nique \2t2\ . Roughly speaking, the elliptic reconstruction technique, as far as energy 
estimates are concerned, allows to neatly separate the time discretization analysis 
form the spatial one. This technique, which has been adapted to tackle fully-discrete 
schemes via energy methods for conforming methods [29], is extended in this work 
to the nonconforming setting, to all methods that meet the two requirements above. 



Briefly said, the idea of elliptic reconstruction — denoting by u the solution of (1.3) 
and by U that of the discrete problem — consists in building an auxiliary function 
w, called the elliptic reconstruction of U , which satisfies two key properties: (a) a 
PDE-like relation binds the parabolic error u — w with data quantities only involv- 
ing w — U and the problem's data, /, a, and Uq, (b) the function U is the Ritz 
projection of w onto Sh- Note that w is an analysis-only device that, despite its 
name, it is not a computable object. Fortunately, computing w is not needed in 
practice, as it does not appear in the resulting a posteriori bounds. 

We believe that it is possible to obtain similar a posteriori error estimates, for 
each single method at hand, by working directly, i.e., without using an elliptic re- 
construction technique, but this will inevitably lead to further complications which 
may render the analysis quite involved, especially for the fully discrete scheme. This 
prejudice of ours is testified by the somewhat surprising lack of previous rigorous 
results in the literature. Finally, we point out that it is possible to follow a simi- 
lar approach to ours in order to derive a posteriori error estimates in lower order 
functional spaces such as Loo(0,T; L2(/2)). 

We remark that new estimators arise in the derivation of fully discrete a posteriori 
error bounds, due to the time-dependent diffusion tensor considered in this work, 
compared to |29l where only time-independent diffusion coefficients are addressed. 
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The following is an outline of this article. After introducing the notation and the 
method in ^ the elliptic reconstruction is used to develop an abstract framework 
for spatially semidiscrete schemes in Sj3] and their fully discrete counterpart in ^ 
The actual error estimators for each particular method are then consequences of our 
abstract framework and specific elliptic error estimators such as the ones presented 
in [21 I2S1 HH mi 123 [H] • Moreover, i n jpj we prove a posteriori bounds for the cor- 
responding steady state problem of ( |1.1| ) for IPDG, thus extending existing results 
[H [26l [231 ES] to the case of general (non-diagonal) diffusion tensor, with minimal 
regularity assumptions on the exact solution [THl cf.]. These a posteriori bounds 
are then combined with the general framework presented in in Sj3l and fS] to deduce 
fully computable bounds for the DG-approximation error of the parabolic problem. 
Last in ^ we summarize results from computer experiments aimed at exhibiting 
the reliability (derived theoretically) and efficiency of the error estimators in the 
special case of the IPDG method. 

2. Preliminaries 

2.1. Functional analysis tools. Given an open subset lo C R'', we denote by 
U'{(jj), 1 < p < oo, the Lebesgue spaces of functions with summable p-powers on 
Lu. The corresponding norms || • ||Lp(tj); the norm of L2(w) will be denoted by ||-||i^ 
for brevity; by (•, •)^ we write the standard L2-inner product on u. When lj ~ fl 
we omit the subindex. 

We denote by H^(aj), the standard Hilbert Sobolev space of index s G R; HQ(a;) 
signifies the subspace of H (w) of functions with vanishing trace on the boundary 
duj. The Poincare-Friedrichs inequality 

\\v\\n < CvY\\Vv\\n for v e Hj(f2) (2.1) 

turns the H^(]7) seminorm into a norm on Hj(f2). We consider thus ||V-||i7 to be 
the norm on Hj(i7). We will use also H~^(i7), the dual space of Hg(J7), equipped 
with the duality brackets (• | •). Namely, if / G H~^(i7) then for each </> G HJ(J7) its 
value on <j) is denoted by {g \ cj)) which coincides with (5,<^) if <? € L2(w). Thus the 
norm of g is given by 

hU-^iO)--- sup ll^. (2.2) 

0eHi(r2)^{o} l|V0|| 

The duality pairing (H^^,HJ) allows us to define, for each fixed t G (0,r] the 
elliptic operator s/{t) : HJ(J?) — > H^^(i7) where 

(j^(i)w I (/)) := {a{t)\/v, \7(j)) (^ f {a{x, i)Vi;(x)) • \7(j){x) dx\ (2.3) 

for (j),v & Hj(i7). Here, and throughout the paper, we use the shorthand f{t) = 
f{-,t), for a function / : [0, T] x J? — > R. Note that thanks to the uniform parabolic 
assumption given by (2.13[) and the Lax-Milgram Theorem, the definition of the 



operator s^{t) is well defined and yields an isomorphism between H~ [Q) and 
Hg(j7) [21]. In other words, the operator s/lt) induces a bounded coercive bilinear 
form 

(w, w) G Hj(f2) X Hj(i7) ^ {.s/{t)v I u;) G R, (2.4) 

which we will extend later to a larger nonconforming space. We stress from the 



outset that although the bilinear form in (2.4) will be extended later to larger spaces. 



the operator ^(t) will not and it will only act on Hq(!7) functions throughout the 
discussion. 



A POSTERIORI CONTROL IN PARABOLIC DISCONTINUOUS GALERKIN 5 

For 1 < p < +00, we also define the spaces LP{0,T,X), with X being a real 
Banach space with norm || • || x , consisting of all measurable functions v : [0,T] ^ X, 
for which 

\Hlp{o,T:X) -^ [ ll«Wllxdtj < +00, for l<p<+oo, 

||^'||L~(o.T;jf) :==csssupo<t<y ||u(t)||x < +00, for p = +00. 

Finally, we denote by C{0,T;X) the space of continuous functions v : [0,T] — > X 
with norm ||w||c(o,T;X) := maxo<t<T lk(i)||x < +00. 

2.2. Finite element spaces. Let 5^ be a subdivision of H into disjoint open sets, 
which we call elements. We assume ^ to be parametrized by mappings F^, for 
each K € £^, where i^^ : k — > k is a diffeomorphism and k is the reference element 
or reference square. The above mappings are such that /? — Uk^stk. We often 
use the word mesh for subdivision, and we say that a mesh is regular if it has no 
hanging nodes; otherwise the mesh is irregular. Unless otherwise stated, we allow 
the mesh to be 1-irregular, i.e., for d ~ 2, there is at most one hanging node per 
edge, typically its center; for d = 3 a corresponding concept is available. 

For an integer p > 1, we denote by ^p{k), the set of all polynomials on k of 
degree p, if k is the reference simplex, or of degree p in each coordinate direction, 
if k is the reference cube. We consider the discontinuous Galerkin finite element 
space 

S = SP{£r) := {v e L2(r2) -.v^oF^e ^p{k), n e .^}. (2.6) 

By F we denote the union of all sides of the elements of the subdivision ^ 
(including the boundary sides). We think of F as the union of two disjoint subsets 
F = Fq\J /Int, where Fq is the union of all boundary sides. 

Let two elements i<^ ^ kT G -S^ have a common a side e :— R'^ Ok^ C /Int- Define 
the outward normal unit vectors n+ and n on e corresponding to di<^ and 9k~ , 
respectively. For functions g : i? — > IR and <p : Q ^ V^ uniformly continuous on 
each of k , but possibly discontinuous across e, we define the following quantities. 
For q+ := q\di^+, q^ := gja^- and ^+ := (j)\oi^+, (f>^ ■= ^Iok-, we set 

We :- liq+ + q-), We := ^(0+ + c/-"), ^^ 7) 

M^ := q+n+ + q-n-, 101^ := 0+ • n+ + 0- • n"; 

if e is a boundary side (e C Fq) these definitions are modified to 

{g}e :- q+, me ■■= C/.+ , M, := q+n+, ^^ := cj)+ ■ n+ . (2.8) 

We introduce the mesh-size as the function /i : i7 — > R, by h(x) — diaruK, if x G k 
and h{x) = {h}, if a; G -T. The shape-regularity of the subdivision £^ is defined as 

fi{,^) := sup ^, (2.9) 

where r^ is the radius of the largest ball that fits entirely in k. 

We shall use the gradient's regular part operator, V^, of an elementwise differ- 
entiable function v defined by 

(V^f)U:=V(«U) forKG ^. (2.10) 

Note that the full distributional gradient, Vw, consists of an extra term taking into 
account the jumps of v across the edges (with a — sign for historic reasons): 

Wv = y^v-M5r,^,, (2.11) 
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with dpi^^ denoting the Dirac distribution on the interior skeleton /Int- Finally, we 
consider some shorthand notation for quantities involving the diffusion tensor a. In 
particular, we define the elenientwise constant functions at,, aj : ^ x [0,T] — >■ IR by 

a»(-,i)U:=IIIVa(-,i)|2|lioo(,) and ai,i;t)\^ := \\\iVai;t))-%\\li^^^, (2.12) 

for K e £^, and aj — {aj}, a\, =^ ({l/a[,})^^, on F, where | • I2 denotes the Euclidean- 
induced matrix norm. Finally, let 

ati(<) := maxaj(a;,i) and a\,{t) — mm a\,{x,t). (2-13) 

2.3. Spatial DG discretization. Introduce the DG space 5^ := 5 + Hq(]7), and 
a corresponding DG bilinear form B : S^ x ,y ^ R, which we assume to be an 



extension of the bilinear form defined by (2.4 1, viz 



B{t; V, z) = {s/{t)v I z) for all v,ze HJ(I2), t e (0, T] . (2.14) 

Though B is time-dependent, we do not write it explicitly in the semidiscrete case 
and omit the t. 

The space S^ is equipped with a, DG norm, denoted |||'||| and depending on the 
method at hand, which extends the energy norm, i.e.. 



|||t,||| ^ \\y/a{-,t)Wv\\ ioi- alive Rl(f2), (2.15) 

for t € [0, T]. Also here, the norm is time-dependent, but this dependence is not ex- 
plicitly written. A norm equivalence between the energy norm |||'||| and || •\/a(-, i)V-|| 
in Hj(i7), uniformly with respect to t, suffices for all the bounds presented below 
to hold, modulo a multiplicative constant; but, we eschew this much generality for 
clarity's sake. 



The semidiscrete DG method in space for problem (1.3 1, reads as follows: 

Find U e C"'i(0,T; Sh) such that 

(2.16) 

{dtU,V) + B{U,V) ^ {f,V) iorVeS,te[0,T]. 

We stress that these assumptions are satisfied by many DG methods for second 
order elliptic problems available in the literature, possibly by using inconsistent 



formulations |3]. Moreover, (2.14) and (2.15) are satisfied by IPDG (along with the 



corresponding e nergy norm) considered below as a paradigm. 

Assumption (2.141 implies the consistency of the bilinear form B on IlJ(i?), i.e. 



{dtU,v)+B{u,v)^ {f,v), for aU w e Hj(i7), (2.17) 



where u is the exact (weak) solution to the initial-boundary value problem ( |l.3[ ). 
Note that this is the same as writing 

dtu + s^u^f. (2.18) 

2.4. Fully discrete solution. To further discretize in time, consider an increasing 
time partition {in}„=o ^v' ^^'^ ^^^ corresponding time-steps t„ = t„ — i„_i, for 
n = 1, . . . , TV. For each n = 0, . . . , A^, we consider that S" is a DG finite element 
space of fixed degree p built on a partition ^n, which may be different from ^n-i 



when n > 1. In [5.3 we will say more about the sequence of meshes and the 
compatibility relations among them. 

Let /"(x) := f{x,tn) and let f/" be the projection (or an interpolation) of uq onto 
the finite element space 5°. We say that {C^"}„^o at is a fully discrete solution of 



( 1.3 ) if, for each n = 1, . . . , A^ we have that [/" S iS" satisfies 

((C/"-t/"-i)/T„,T/)+S"(C/",V^) = (/",t/) for aU ye 5", (2.19) 
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Since the elliptic operator si (and the bilinear form B) depend on time, in the 
fully discrete setting, we denote their value at time t by j2/(t) and B(£), respectively 
and when t = tn we take jz/" := ^(in) and i?" := B{tn). 

Noting that the term [/""^ can be replaced by 7T"C7"-\ where 77" : L2(/?) -)■ 5" 
is the orthogonal projection, we consider a slightly more general situation where 



jjnjjn-i jj^ ( |2.i9[ ) is replaced by /"[/""i; here /" : S"""! ^ S"" is a general date 
transfer operator, depending on the particular implementation. The operator J" 
may coincide with 77", but it may be an interpolation operator for example. The 
general fully discrete Euler scheme then reads 

(([/"- 7"C/"-i)/T„,y)+B"(f/",y) = (/",V^) for all ye 5". (2.20) 

We have taken /" = /(tn), but we could take a more general approximation than 
f{tn), for example, a good choice is also given by /" := J^ " /(s) ds, for which a 
suitable modification of our arguments leads to similar results. 

3. Abstract a posteriori bounds for the semidiscrete problem 
We derive next an abstract a posteriori error bound for the quantity 

\\u-Uh,^o,T;^y.= (^J \\\uit,-)-U{t,-)\ry^\ (3.1) 

where |||'||| denotes the appropriate (space) energy norm. 

In the a posteriori error analysis below, we shall make use of the idea of ellip- 
tic reconstruction operators introduced in [32 for the semidiscrete problem and 
extended to fully discrete (conforming-in-space) methods in ^\. 

3.1. Definition (elliptic reconstruction and discrete operator). Let U be the 



(semidiscrete) DG solution to the problem (2.161. We define the elliptic reconstruc- 
tion w e HQ(f2) of [/ to be the solution of the elliptic problem 

B{w,v) = {AU -nf + f,v) for ah w e HJ(J7), (3.2) 

where 77 : L2(J7) -^ S denotes the orthogonal L2-projection on the finite element 
space S, and A : S ^ S denotes the discrete DG operator defined by 

{AZ, V) = B{Z, V) for all V e S, (3.3) 

for each Z G S. Note that this is valid on [0, T]. 

3.2. Remark (the role of the elliptic reconstruction). The elliptic reconstruction 
is well defined. Indeed, AU d S is the unique L2-Riesz representation of a linear 
functional on the finite dimensional space S^ and the existence and uniqueness 



of (weak) solution of (3.2), with data AU — 11 f + / G L2(i7), follows from the 



Lax-Milgram Theorem in view of (2.141 



The key property of w is that the DG solution U of the semidiscrete time- 



dependent problem (2.16) is also the DG solution of the steady-state boundary- 
value problem (3.2). Indeed, let W G S be the DG-approximation to w, defined by 
the finite dimensional linear system 

B{W,V)^{AU-nf + f,V), (3.4) 

for ah V e S, which implies B{W,V) = {AU,V) = B{U,V) for all V e S, i.e., 

W :^U. 

3.3. Definition (error, elliptic and parabolic parts). We shall decompose the 
error as follows: 

e := U ~ u — p — e, where e := w — U, and p :— w — u, (3-5) 

where w — w{t) denotes the elliptic reconstruction oi U = U{t) at time t E [0,T]. 
We call e the elliptic error and p the parabolic error. 
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3.4. Lemma (semidiscrete error relation). Let u be the solution of Problcm \l . 1\ 
U denote the solution of the DG scheme {2.16). Then, we have 



{dte, v) + B{p, v)^0 for all v e HJ(I2). 

Proof. For each v G HJ(J?) we have 

{dte, v) + B{p, v) ^{dtU, v) + B{w, v) - (/, v) 

= {dtU, v) + {AU -nf + f,v)- if, v) 

={dtU, nv) + {AU, nv) - (/, nv) = o, 



(3.6) 



(3.7) 



where in the first equality we used (2.17); in the second and fourth equalities, we 



made use of Definition |3.1[ in the third equality the properties of the orthogonal 

D 



L2-projection onto S are used; finally, the last equality follows from (2.161. 



3.5. Definition (conforming-nonconforming decomposition). In the theory de- 
veloped below, we shall consider the decomposition of the DG solution U & S into 
conforming (continuous) and nonconforming (discontinuous) parts as follows 

U^Uc + Ud, (3.8) 

where Uc € Sc := Hj(f2) n S and Ud := U - Uc e S. Note that at this point 
we do not specify any particular decomposition, thus keeping the choice of such a 
decomposition at our disposal. Let 

e^ :=e-t/rf = [/c-weHj(J7), and e^ := e + C/rf = w - C/^ G Hj(f2). (3.9) 



3.6. Theorem (long-time a posteriori error bound for DG). Let u and U he the 
exact weak solution of ( 1.3) and the DG solution of the problem ^2.16 ), respectively. 
Let w be the elliptic reconstruction ofU as in Definition\3.1\ Assuming {2.14) and 



{2.15), and that a decomposition of the form {3.8) is available, then the following 
error bound holds 

li^-u|lL2(0,T;,r) <3|1w-C/||l2(0,T;,5-) +2||C/d||L2(0,T;^) 

+ 2\\dtUd/V^^L,^o,T-M-\n))+nuo-U,{0)l 



Proof. Set V = Be in (3.6); then, in view of ( |3.9[ ), we have 

{dte,, e,) + B{p, p) - Bip, e,) - {dtUd, e,). (3.11) 

Recalling (2.14), (2.15) and that p,ec e HJ(J7) for every t S [0,r], using the 
Cauchy-Schwarz inequality, and the duality pairing (H^ , Hq) , we arrive to 

1 



:djee||^ + |||p|||^<|||p||| 



Also, (3.91 implies 



Hi(r2) 



< 



\\dtUd\\n-^[n)\ 

IIpIID/v^- 



VeJ 



(3.12) 



(3.13) 



Setting /i 



which implies 



fc 



:~ \\dtUd/ ^/oi^\\Yi-i(Q) in (3.12) and rearranging lead to 



IHII 



< 



\pmi 



h 



hh, 



dt\\e,r+\\\p\r<Aii!+ii). 

Integration on [0, T] and taking square roots yields 

IIpIIl2(o,T;j^) < l|ec(0)|| +'2\\ec\\L2{o.T;.y)+'2\\dtUd/V^\\L2{o,T;H-^f2))- 



The assertion follows using triangle inequality on (3.51 and (3.9 1. 



(3.14) 
(3.15) 

(3.16) 

D 
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3.7. Theorem (short-time a posteriori error bound for DG). Let the assump- 
tions of Theorem |3.6| hold. Then the following error bound holds: 



(3.17) 



\\U - U\\L2{0,T;y) <(V3/2 + 1)||W - C/||l2(0,T;J^) + V^/'^\\Ud\\u(0,T;J^) 

+ ndtUd\\LHo,T-Mif2))+V2\\uo-UMl 
Proof. Let Tq e [0, T] be such that 

||e,(To)||= max^||e,(t)||=:ii;e. 



Then, (3.11) impUes 



2dt|lec|l^ 



Illpr<||lp||l niacin +i?c|15,C/,|l 



which, after integration on [0,ro], yields 



:E' 



or 



II/'IIl2(0,To;J?') <2ll'2c(0)|| + ||p||L2(0,To;J5*)l|ec||L2(0,To;J?') 
+ Ec\\dtUd\\L^0,To;L2(n)), 

1, 



V< lHWnM,2 



|ec(0)|| 



'^c\\L2(0,T-,y) + ll^tf^d||i,i(0,T;L2(i7))- 



(3.18) 



(3.19) 



(3.20) 



(3.21) 



Going back to (3.19), upon integration with respect to t between [0,r], we obtain 
1, 



2.IpIIl2(o,T;^) < ^l|ec(0)|P + ^||ec||L2(o,T;^) + 4^' 



\dtU4^ 



Li(0,T;L2(r2))' 



(3.22) 



which, in conjunction with (3.21) gives 



IIpIIl2(o,T;^) < 2||ec(0)||' + ^||ec|lL(o,T;^) + M\dtUd\\h^o,T-Mmy (3-23) 



the final bound now follows using the triangle inequality on (3.5 1 and (3.9) 



D 



3.8. Remark (long- versu s sho rt-tim e bounds). We note that the crucial differ- 
ence between bounds (3.10) and (3.17) is that in the latter the L-'^-accumulation 
term \\dtUd\\]^i(^Q^T;,s^) is present; this implies 

\\dtUd\\LHo,T-M(m < ^^l|9tf/<i||L2(o,T;L2(r2)), (3.24) 

which may be preferable if T < 1, but can be inefficient for long-time integration. 
On the other hand, the corresponding term in (3.10) is ||9tC/d/-v/Q^llL2(o,T;H-i(^))' 



which can be a bit less inefficient when the diffusion tensor a varies substantially 
on f2. We note, however, that it is possible to avoid dividing by the factor a\,, 
by equipping H^ (J7) with the dual norm of the energy norm (2.15) in Hq(^). 



In practice, however, this improvement is relevant only if the dual norm is calcu- 
lated explicitly [30^ . Alternatively, one can apply a Poincare-Friedrichs inequality 
to bound the dual norm by the L2-norm, which results into the reappearance of the 
factor at. 



3.9. Remark (elliptic a posteriori error estimates). The bounds (3.10) and (3.17) 



are not (yet) explicitly a posteriori bounds: ||w — {/||l2(o,t;^) still needs to be 
bounded by a computable quantity. To this end, given g S L2, consider the el- 
liptic problem: 

find z £ Hq(J7) such that 
— V • (aVz) ~ g in [2, 2; = on Oil, 



(3.25) 
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whose solution can be approximated by the following DG method; 

find Z (z S such that 

B{Z, V) = (5, V) for ah V £ S. 

If assume that an a posteriori estimator functional S exists, i.e., 

\\\z-Z\\\ <S(Z,a,g,3r), 



(3.26) 



(3.27) 



then we can computably bound ||w — U\\i^^^qx-j^) in (3.10) and (3.17) through 



lk-C^llL2(o,T;,y) 



< 



f S{U,a,AU -nf + f,.!7)'' 
Jo 



1/2 



(3.28) 



A posteriori bounds for various DG methods have been studied, under different 
assumptions on data and admissible finite element spaces, by many authors |6l [26l 



ESI El EH mi [ig. Thus Theorems [S^e] and [3J] can be app lied to any DG— and 
more generally to any non-conforming — method satisfying (2.14) and (2.15), and 
for which (3.27) is available. The object of i4]is to address this for IPDG. 



4. A POSTERIORI ERROR BOUNDS FOR THE INTERIOR PENALTY DG METHOD 

Here we extend the energy-norm a posteriori bounds for the family IPDG meth- 
ods cf. 13 EHl El] for the Poisson problem, to the case of the general diffusion prob- 
lem ( |3.25| ). A similar analysis has recently appeared also in [19 , while a related 
DG method based on weighted averages for anisotropic and high-contrast diffusion 
problems can be found in 20J. We stress that our results can be generalized as to 
allow for inhomogeneous or mixed boundary conditions following |23ll26l resp.]. 

4.1. Definition (IPDG method). Fo r z,v & y, the bilinear form B : ,yx,y ^R 
for the IPDG method for the problem (3.25) can be written as 

B{z,v):^ f {aW^z)-V^v+ f (^{ailVt;}- H-{a7TVz}- H+^r H -H ), (4.1) 
Jn J r 

for 9 e {-1,0,1}, where 77 : [L2(i7)]'' -^ S"* denotes also the orthogonal L2- 
projection operator onto S'^ , and the penalty function cr : 7^ —> R is defined by 



C„ 



.m(^) 



{aj} 



where the constant C„ 



,A'(^) 



(4.2) 



> depends on the shape-regularity of the mesh i^ 



and on the smallest possible Cq > 1 such that 



c- 



< 






<a, 



(4.3) 



for every pair of elements k^ and k sharing a common side. 

We also define the corresponding energy norm |||'||| for the IPDG method by 



M 



(wV^y^vf + WV^Hrr) 



1/2 



(4.4) 



for V £ S^ . Note that for w e S*, we have 77Vw — Vw and, therefore, B can be 
reduced to the more familiar form 

73(z, v) := / [a^^z) ■V^v+ f {9{aVv} ■ {zj - {aVz} • H + ^ W ' H ) , (4-5) 
Jn Jr 

for z, w e 5 [SmilMl El cf.]. Observe that both ( pUj ) and ( [2l5| ) hold for the 
particular B and |||-||| defined above. 
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4.2. Remark (conforming part of a nonconforming finite element function). The 
space-discontinuous finite element space S contains the conforming (continuous) fi- 
nite element space Sc = Sr\liQ{f2) as a subspace. The approximation of functions in 
S by functions in Sc will play an important role in our derivation of the a posteriori 
bounds. This can be quantified in the following result, which is an extension of [27l 
Thm. 2.1]. For other similar results we refer to [iOl [51 E51 [TU] . 

4.3. Lemma (bounding the nonconforming part via jumps). Suppose ^ is a 
regular mesh and a is elementwise (weakly) difFerentiahle. Then, for any function 
Z G S there exists a function Z^ G Sc such that 

\\Z-Zc\\<C4VhlZ}\\r, (4.6) 

and 

\\V^V^{Z-Zc)\\<C2\\V^m\\r, (4.7) 

where Ci , C2 > constants depending on the shape-regularity, on the maximum 
polynomial degree of the local basis and on Ca- 

The proof, omitted here, follows closely that of f2(r, Thm. 2.2]. Lemma 4.3 can 
be proved for irregular (i.e., with hanging-nodes) meshes [26( Thm. 2.3], in which 
case Ci and C2 depend on the maximum refinement and coarsening levels imax- 

4.4. Lemma (a posteriori bounds for IPDG method for elliptic problem). Let 
=5^ be a regular and a is elementwise (weakly) differentiable. Let z and Z be given 



by j3^ and ^3J§. Then 

\\\z-Z\\\<^ip(Z,a,g,^), (4.8) 

where 

^ip{Z,a,g, 3r) -.^ CKf (\\h/y^,{g + V,^ ■ {aV^Zm 

(4.9) 

|v/Va^IaV^^lllr.„, + IIA 

and K'^ := maxf^ \fo^/oA,, where C > depends only on ni^) and Ca- 

Proof. Our proof is inspired by [261 123| . Denoting by Zc G Sc the conforming part 
of W as in Lemma 1431 we have 



e := z — Z ^ Cc + Cd, where Bc := z — Zc, and e^ := Zc — Z, (4-10) 

yielding Cc € HJ(J?). Thus, we have B(z,ec) = {g,ec)- Let TIq : ^2{^) — >■ R denote 
the orthogonal L2-projection onto the elementwise constant functions; then Tloec € 
S and we define r] := Cc — Iloec. 
We also have 

B(e, ec) = B(z, Cc) - B(Z, ej = [g, e,) - B{Z, 77) - B{Z, TToe,) = {g, ri) - B{Z, 77), 

(4.11) 
which implies 

ll^VedP =B(e„e,) = {g,ij) - B(Z,rj) - B{ed,ec). (4.12) 

For the last term on the right-hand side of ( |4.12[ ), we have 

\B{ed,ec)\ < \\y/^V3red\\\WaVec\\+-J2 Yl at)UII^(^Ve,)U||e|| M /^ll^, 



2 



(4.13) 
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where k+ and k^ are the (generic) elements having e as common side. Using the 
inverse estimate of the form ||-\//i^||e < C'||V^||k for V = II\7ec, and the stabihty of 
the L2-projection, we arrive to 

|S(erf,e,)| < WVaV^edWWV^Ve.W + CK^\\VaVe,\\\\V^le4\\r. (4.14) 



Finally, noting that Jcd] = |^] , and making use of ( |4.7| ) we conclude that 

\B{ed,e,)\ < CKfWVaVecW \\V^lW}\\r. (4.15) 



To bound the first two terms on the right-hand side of (4.12), we begin by an 
elementwise integration by parts yielding 



(g, 7/) - BiZ, ry) = / (5 + V^ • (aV^Z))r/ - / {ry} laV Zj ds 

+ J e{anVrj} -iZjds- J a iZj ' M ds- 
The first term on the right-hand side of (4.16) can be bounded as follows: 

/ (g + V • (aVZ)) 7/ < \\h/^^^{g + V^ ■ iaV^Z))\\\\y^h-^Tj\\; 
upon observing that ||/i^-'^?;||k < Cj|Vec||K, this becomes 

/ (.g + V • (aVZ)) 77 < CK-^\\h/ ^,{g + V:^ ■ {aV^Z))\\\\^Ve,\\. 



(4.16) 



(4.17) 



(4.18) 



For the second term on the right-hand side of ( 4.16[ ), we use a trace estimate, the 
bound ||/i~^7/||„ < CIIVecIIk and we observe that Vry = Vgc, to deduce 



MlaVZlds < C/d|^/SVe,|| \\^h/^aay^Z}\\r,^ 



-Tin 



(4.19) 



For the third term on the right-hand side of ( |4.16 ), we use Vjy — Ve^ and, 
working alike to (|4.13[), we obtain 



6l{ai7Vry} 



<CK-^\6\\\^Ve,\\\\^ 



and finally, for the last term on the right-hand side of (4.16), we get 



<CA'f||V^V£ 



The result follows combining the above relations. 



^im 



(4.20) 



(4.21) 



D 



4.5. Theorem (a posteriori bounds for IPDG method for parabolic problem). 
Let u, U be the exact weak solution of { 1.3), and the IPDG solution of the problem 
{2.16), respectively, and let a be elementwise (weakly) diffcrcntiablc. Then, the 
following error bound holds: 

( 0I (TJ ^ ATT TT £ I £ 1^ \ i „— 111. /TT ^ '. 



Ik - C^llL(o,T;^) <G f (<^il((7, a, AU-nf + /, ^) + a;'\\VhldtUjfr) 

JQ 

-fC(||wo-C/(0)||+||\/^l 



(4.22) 
If we assume also that u,U e C(0, T; HJ(I2)) n H^(0, T; L2(f2)), then the following 
bound also holds: 

\W-U\\l,io,T;^)<C j Sl{U,a,AU^nf + f,:7) + ci^j \\^/hldtUl\\rf 



c(||uo-c/(o)|| + 11 a 



(4.23) 
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Lemma 4.4 in conjunction with (4.61 



Proof. The results follow immediately from combining Theorems 3.6 and 3.7 with 

D 



Finally, we give a result on useful properties of the IPDG bilinear form and of 
the norm, which will be useful in ^ 

4.6. Lemma (continuity of B and stability of the L2-projection). Consider the 
notation of ^ and let B and III 'III denote the IPDG bilinear form 14.51) and the 



DG-norm {4.4). Then for Z,V £ Sh we have 

BiZ,V)<CK^\\\Z\\\\\\V\\\. 
Moreover, for v G Hj!,(i7), the L2 -projection is DG-norm-stable, i.e., 



\nv\\\ < CK. 



sr\ 



(4.24) 



(4.25) 



Proof. We omit the proof of (4.24) which mimics that of (4.15). For stability, note 

\\nv\f ^WaV^iBv - nov)f + \\^lv- 7Tt;]||f, 

<c{\\alh-^{nv - nQv)f + \\a\h-'^{v ~ nv)f + WafVsriv - nv)f\ 



<C[\\alh-\v - nov)\\' + ||aiV^(t 
<C||(as/aO^\/aV,Tv|P, 



which implies (4.25). 



(4.26) 

D 



5. A POSTERIORI ERROR BOUND FOR THE FULLY DISCRETE SCHEME 

In this section we discuss the abstract error analysis for the fully discrete scheme 



defined in ^2.4 



5.1. The elliptic reconstruction and the basic error relation. Extend the 
sequence {f/"} into a continuous piecewise linear function of time: 

C/(0) = [/° and U{t) = /„(t)C/" + /„_i(t)t/"-\ (5.1) 

for t € [in-ijin], and n = 1, . . . ,iV, where the functions In and In-i are the La- 
grange basis functions 



ln{t) .= l[t„_i,t„] + — l[t„,t„+i 



(5.2) 



Using these time extensions and the (time and mesh dependent) discrete elliptic 



(5.3) 
(5.4) 



operator A" of Definition |3.1| with respect to S*" , defined by 

A"Z e 5" such that (A"Z, V) = B'\Z, V) for all V G S"' 



we can write the scheme (2.20) in the following pointwise form: 
dtU{t) + A"U" == (/"[/"-I - U"-^)/Tn + n'^f", 

for all t e {tn-i,tn) ,n — 1, . . . ,N. 

We like to warn at this point that we use the same symbol U{t) to indicate the 
fully discrete solution time- extension in this section, and the semidiscrete solution 
in SQ This should cause no confusion as long as the two cases are kept in separate 
sections. 

For each fixed t G [0,r], and the corresponding n = 1,...,N such that t G 
(t„_i,i„], we define the time- dependent elliptic reconstruction to be the function 
w{t) e Hj(i7) satisfying 

w{t) = ln{t)w" + ln-i{t)wl-\ (5.5) 
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where w" G Hq(^) is the elliptic reconstruction of C/" defined impUcitly as the 
(weak) solution of the elhptic problem with data A^U"', i.e., w" satisfies 

^"u;"=A"C/", (5.6) 

and w"^ is the forward elliptic reconstruction of C/"^^, defined as the solution of 
the problem 

.a/^'-^wl-^ = A^-^r'V^-^, (5.7) 

where the operator A"^^ : S*" — > 5*" is defined by 

A^-^Z e 5*" such that (A^-^Z, V) = Bl-\Z, V) for aU V e S'\ (5.8) 

i?"^ being the nonconforming bilinear form corresponding to ^"~^, but with 
respect to the space 5" (in contrast to i?"^^ which is defined with respect to 
5"^^). For instance, for IPDG, we have 

Bl-\Z,V):= J2 f{a{tn-i)^Z)-^V 

{e{a{tn-i)nnWV} ■ IZj - {a(t„_i)77„VZ} • 1^1 + a„ IZj -IVj). (5.9) 

Using this definition of w on [i„_i,i„], the equation ( |5.4[ ), implies 

dtUit) + £/{t)wit) = {rU''-^ - C/""^)/t„ + /" + s^it)wit) ~ ^"u.". (5.10) 

Subtracting the exact equation from this identity we obtain 

dt[U -u]+ £/[w -u]^ (/"[/"-! - f/"-i)/r„ +(/" - /) +{£/w - i/"u;") (5.11) 

for all t G {tn-i,tn) and n — 1, . . . ,N, this leads to the following technical basis of 
this section. 

5.2. Lemma (fully discrete error relation). With the notation introduced in this 
section, let e = U — u (fuU error), p := w — u (parabolic error) and e := w — U 
(elliptic error). Then we have 

dtc + £/p^ {rU''-'^ ~ U''-^)/Tn + /" - f + s^w- ^"u.", (5.12) 

on [tn-i, tn], for all n = 1, . . . , N. 



L 



Proof. Replace the new notation for the errors in (5.11). _ 



5.3. Mesh interaction, DG spaces and decomposition. The domain ^'s sub- 
divisions (also known as meshes) {^n}n=i3 n ^''^ assumed to be compatible in 
the sense that for any two consecutive meshes, say ^ and ^n-i, we have that 
c5^ is a constructed from S^n-i in two main steps: (1) =5^_i is locally coarsened by 
merging a chosen subset of elements then (2) the resulting coarsened mesh is locally 
refined [291 [30]. This procedure leads to meshes which are locally a refinement of one 
another. For example, in the following diagram the mesh S/'n-i bas some elements 
marked (in red) for coarsening and (in blue) for refinement and the mesh ^ can 
be thus obtained in the two steps: 



refine 



or _ LAJaiocii acinic y_ _ ^ 

LJ±^ ^U ^U (5.13) 

For each n — I, . . . , N, we denote by ^n the coarsest common refinement of 3^„_i 
and 3^n- The finite element space corresponding to =3^ being 5", we shall be using 
the space Sn which is the finite element space with respect to ^n- Furthermore 
we denote by ^„ :— 5" + Hq(J7) and by ,y := X]n=o "^"' ^^^ minimal space that 
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contains all these spaces. These spaces are equipped with the same type of norms 
given as in Section (23) 

The conforming-nonconforming decomposition of U, that we shall be using is 
performed as follows: 

(a) For each given i„, with n=l....,A^— Iwe assume that the following two 
decompositions exist for [/", 

[/" = L/? + [/" with respect to the mesh ^, 

! (5.14) 

and [/" = UJl^ + C/"_^ with respect to the mesh ■%,+!. 

Note that if the mesh changes, in general, U^ and U^, need not be equal functions. 

(b) For each t S (in_i, in], n — 1, . . . ,N, we define 

Ud{t) := ln-i{t)U^;' + ln{t)U2 and U,{t) :== ln-i(t)Uc+^ + 'n(i)C/"- (5-15) 

5.4. Definition (a posteriori error indicators). We set here some notation that is 
useful to state the main results concisely. We make some assumptions in the process. 
For each time interval [tn-i, ^n], with n = 1, . . . , N, we introduce a posteriori error 
indicators as follows. 

(a) We assume that there exist Cds, C'dgc > such that 

lli7%L<Ceis|klL, for alU- e Hj(l?), (5.16) 

and 

B^{Z,V)<Cd,,\\Z\\JV\l, foranZ,Fe5". (5.17) 

The time-stepping indicator is given by 

^^^_Celsaigc||^„^„_l_^„||^^ (5.18) 

(b) The time data-approxim,ation indicator is 

(c) The mesh-change (or coarsening) indicator is defined as 

1.:=^^ "-JhlZim 1 1] . (5.20) 

(d) The parabolic nonconforming part indicator is given by 

1/2 






„,.._ 21™ ir 1 

and the elliptic nonconforming part indicator defined as 

- ,Xl/2 



(5.21) 



Sn:=[\\Un:+ K;' J . (5.22) 

(e) The space (or elliptic) error indicator is given by 

e„ :=^([/",a(i„),A"C/",5;), (5.23) 

where (^ is a particular choice of an energy-norm elliptic error estimator for the 

given spatial method. Furthermore the forward elliptic error indicator, due to mesh 
change, is given by 

4-1 := ^(/"C/"-\a(i„-i),Ari/"C/"-i,^0- (5.24) 
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(f ) Consider first the auxiliary function of time 



^a,n[S) : — 



ya{s)a{t„) ^ - ya{s) a(t„) 



se[t„„i,t„] (5.25) 



L^(17) 



wliere tfie inner matrix norm is the Euclidean-induced one. This definition is possi- 
ble thanks to a's being symmetric positive definite. The function Xa,n is identically 
zero if the operator is time-independent, otherwise it acts like the numerator of 
a's normalized Holder-continuity ratio.) Then we may define the following operator 
approximation indicators 



c--= 



IJ^njjr 



lH-i(r2) 



1 



% 



1/2 



"n' a.n 



'^ A^ 



Sn ■ 



Mn-lj«[/„-l| 



H-i(r2) 



Ok 



/2 ^2 



t,.-l 



Q := \\[A\-^ - A"] /"[/" 



1 



Ih-1(J7) 






1/2 



1/2 



ab 



(5.26) 

(5.27) 

(5.28) 
(5.29) 



Cn ■— 4n + Cn + Cn • 

(g) Finally, the parabolic nonconforming part indicator of higher order is given by 






(5.30) 



5.5. Remark (computing the H^^(^) norm). The norm H^^(J7) appearing in 
the indicators is easily computable at the cost of inverting a stiffness matrix [301 . 
For many practical purposes, though this has to be replaced by the L2(i7) norm 
times the Poincare-Friedrichs constant Cpp defined in \2.1\ which implies the dual 
inequality 

||w|lH-i(.f2) ^ C'pfII^'II, for all u e L2(X2). (5.31) 

Note that this will not deteriorate most of the indicators. The only indicators that 
may be affected by this change are 5„ and 5„, and it may be possible to provide a 
sharp bound for negative Sobolev norms of [7^, but this seems to remain an open 
question at the time of writing. 

5.6. Remark (computing A^U"^ and similar terms). The operators A" and ^" 
appearing in Definition |5.4[ can be realized in two ways in practice: 



(a) To save time, one can use the fully discrete scheme in pointwise form (5.4 1 to 
evaluate some of these terms. For example 

A"C/" = 77"/" - ([/" - /"C/"-1)/t„. (5.32) 

(b) The corresponding stiffness matrix could be computed and applied to the 
argument. This seems to be necessary for A" . 

5.7. Remark (mesh-change prediction). The mesh-change indicator 7„ can be 
precomputed in a given computation. Indeed, this term does not use explicitly any 



quantity deriving from the solution of the n-th Euler time-step (2.201. This term 



is usually computable when a precise operator J" is available and it involves only 
local matrix-vector operations on each group of elements to be coarsened. 

5.8. Remark (an alternative time-stepping indicator). An equally valid defini- 
tion for the time-stepping estimator 0„ can be given by 
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This alternative definition has the advantage of having no constants, but it is more 
comphcated to compute and it must be reduced to the Li2{f2) norm by using the 
Poincare-Friedrichs inequahty. A good side effect of this ahernative choice is that 
in this case the indicator <^° vanishes; all other estimators remain unchanged. 

5.9. Theorem (abstract a posteriori energy-error bound for Euler-DG). Let 
{C/"}^ be the solution of {2.20) and U its time-extension as defined by {5.1) and w 



the eUiptic reconstruction as defined by (5.5p. Then, with reference to Definition\5.4 
for each m = 1, . . . ,N, we have 



II" - ^llL2(0,t,„;^) < ll"(0) - UciO)\\ + 377p,,„ + V2lJ^^r. 

\ 1/2 



1 



Yl ^"^" 



/ ^ ^n'^rfi 



(5.34) 



where the parabolic-error estimator is defined as 

/ m N 



1/2 



\n=l 



and the elliptic estimator is defined by 



1/2 



(5.35) 



(5.36) 



We spread the proof in paragraphs |5.10| - |5.13l 
5.10. The energy identity. As in the proof of Theorem |3.6| to get an energy 

(5.37) 



identity out of (5.121, we will test with the error's conforming part 

Cc -.^ e - Ud = p + Cc- 



Start with combining (5.12) and definition (5.5) to get 

dtCc + J^p^ dtUd + (/"[/"-I - C/""^)/t„ + /"-/ + ^«; - ^"w" (5.38) 
Testing the above relation with Cc we obtain the following energy identity: 



1 



At\\e£ + \\p\\l={dte,,e,)+B{p,p) 



=B(p,e,) + (9tC/d,e,)+ , 
+ ((/"C/"-i - C/"-1)/t„ + /"-/, ec) 



(5.39) 



Integrating (5.39) from to im G {^^T], for an integer in, 1 < m < N fixed, we 
may write the integral form of the energy identity 



l\\ec{t^)\\' + I "^ \\p\\l^l\\eM\? + 



B{p,ec) 



E 



{£/w - J^^'w" I Cc) 



tn 



{{I'^U" 



[/" 



')/T„ + r-./,e,)) 



(5.40) 



{dtUd,e,) + -Y^{\\u{t^)-Ul\ 



I ^cKj^n 



^0 + ^l(im) + Mtn^) + ^3(im) + Mtm) + M^m)- 
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To obtain the a posteriori error bound for scheme (2.20), we now bound each of 
J^iltm), i = 1, . . . , 5 (^0 needs no bounding) appearing in relation (5.40), in terms 
of either a-posteriori-computable or left-hand-side quantities. 

A term that substantially distinguishes the fully discrete case from the semidis- 
crete one discussed in §3^ is the time-discretization term ^2{tm), so we start by 
bounding this term. 

5.11. Time discretization estimate. To bound J^2(im) we start by working out 
the first factor of the integrand as follows 

Since w" and Cc are both in Hq(J?), we may bound the first term with 



(5.41) 



(K(s) - ^"]zi;" I e,(s)) =. / ((a(s)-a(t„))Vu;")-Ve,(s) 

Jn 

= j (./^)U a{s)a{tr^y^ - ^a{sYa{trM ^a{tn)WwA ■ Ve,(s) 

V^Vee(s) Va(*n)Vu;" 



a{s)a{tn) ^ - Y a(s) \i.{tn) 
Aa,„(s)||^"||J|ec(s)L. 



i^ocin) 



The second factor above can be bounded as follows 

2 



\W 



u;« |u;«) = (v4"C/",ii;") = ||A"C/"||jj-i(^^) ||Vu;"|l 

|U"[/"| 



< 



lH-(^ 



(5.42) 



(5.43) 



where the last step owes to the fact that 

«^liV«;"||'<i?"(u;^«;")-||«;"||^ 
Thus ||w"||^ < ||A"[/"||jj-i(^^)/a^ and we obtain 



Z„(5)([^(s)-^"]w;"|ee(s))ds 



= CnV^I|ec|lL2([t„_i,t„];J^) ■ 

Similarly, we obtain 

" ln-i{s) ([^(s) - ^"-1] u;!^"i I ec(s)) ds 



(5.44) 



1/2 



< 



|^n-lj„f;n-l| 



lH-i(r2) 



t„ 



1/2 



l^cll^ 



(5.45) 



1/2 



ln^l{sf\a,,i-l{sf ds 



t,i 



CnV^I|ec|lL2([t„_i,t„];J^) 



1/2 



(5.46) 
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To estimate the resultant of the integrand's third term in (5.41), we recall the 
elliptic reconstruction's definition and note that in view of (5.4 1 we may write, for 
n > 1 , that 



{s^' 



'-— 1^,,^ — 1 



= ([AIl^-i - A"] /"[/"-\ ec) + (A" [/"C/"-i - [/"] , 7T"ec) (5.47) 



given that A"/"t/""\ A"C/" e S"". The first term on the right-hand side of \hAl\ 
is simply bounded by 



llA\-^ - ^"] /"C/"'\e,) < \{A\-^ - A"] /"C/""iH-i(r2) H^^-l 
and thus, recalling definition ( 5.28[ ), we have 



(5.48) 



(5.49) 



The second term on the right-hand side of (5.47) can be given a simpler expres- 
sion as follows: 



(A"[/"C/"-i - t/"] ,iT"ec) =B"(/"X/"-i - C/",7T"e 



<Cdsc [[/"f/"^ 



[/'^ 



Ii7"e, 



(5.50) 



<CdgcC'cis u"c/" — U"^' 



thanks to the stability of 77" with respect to the energy norm ||-||^ assumed in 
(5.16). Therefore, recalling definition (5.18), we obtain 



Z„_i (A" [/"[/"-I - t/"] , i7"e,) < 0„V^ l|ec|lL,([,„_,.,„];^) (5.51) 

1 

The time error estimate follows: 

•rn 
Mtm) < ^(Cn + 0n) V^(ll PIIl, (t„_, ,t„;^) + II ^c |1l, (t„_„t„;^)) • (5-52) 



5.12. The other error estimates. To bound the spatial error term, ^i(t) in 
(5.40), we simply consider 



^iitm)= f^" B{p,e,)< /"lIpllJlecIL, 



(5.53) 



with the aim of absorbing the first factor in the left-hand side of (5.40) and using 
an elliptic error estimator to bound the second term. 

The term J^3(im) in ( |5.40| which takes into account data approximation: 

^3(t™) = E / ( -^ + /" - /' e- ) • (5-54) 

The first term can be bounded by using the (H^ (J7),Hq(^)) pairing, as we did 
with the time-estimator above: 

E / ((^"C^""' - U^-')/rn + r-f, ee) 

n=l-^«— 1 

m 

< X!(^"+^")^^ll^<=llL2(t„-i,t„;J^)' (^-S^) 



where we have used definitions (5.191 and (5.20). 
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Hence we obtain the bound 

m 
■Mtrn) < ^Wn + 7«) V^(llPll L,(t„_i ,t„;J^) + II ^^ll L, (t„_i ,t„ ;J^)) ■ (5-56) 



We estimate the second-last term on the right-hand side of ( 5.40 1. This term can 
be bounded in two different ways. For concision's sake we expose only the estimate 
that yields smaller accumulation over long integration-times: 



„^i^i„-i 



< Y^ '5«V^(llPllL2(t„„i,t„;J^) + ll<=c|lL2(t„_i,t„; 



(5.57) 



.y) 



by recalling (5.21 1. 



Observing the identity 



\u 



{tn) - U:\\\' - ||ee(i„)|r = \\U2 - U2J' + {U2 ~ C/,"+,e,(t„)) , (5.58) 



we estimate the last term on the right-hand side of (5.401, as fohows: 

^ m—l 

A(U =2 E ( PS - U'S+f + {U2 - t/<?+,ec(t„)) ) 

^O Y\ ('^n'^n + ,^"3-^^ , l|ec(tOII l^nTn) 
Z ^ — ' l<l<7n—l 



(5.59) 



< 



n=l 
m—l 
I ( ^ At,, 



2 1 



- max |lec(t;) 

4 l<l<m~l 



5.13. Concluding the proof of Theorem |5.9[ Combining the energy relation 



( |5.40[ ) with the bounds ( |5.52[ ), ( |5.53| , ( |5.56[ ), ( |5.57[ ) and ( |5.59[ ), we obtain 

-. „ m—l 2 1 

<2 l|e-(0)ll' + i( E '^"^0 + 4 i<?^^^_i ll^^(*')ll' 

m 

+ E(^" + ^" + /3„ + 7„ -I- 5„) v^ l|ec|lL2(t„_i,t„;J?') 

n=l 

771 

+ E((^" + C71 + /3n + 771 + ^n) %/t^ + l|ec|lL2(t„_i,t„;^)) II Pll L^ (0,t„;JJ-) • 



(5.60) 



Choosing m = m^ so that ||ec(tm,)|| = inaxi<;<„j_i ||ec(ti)ll ™ (5.601, yields a 
bound on maxi</<m_i |lec(i;)|| /4, which is then used again to bound the third 
term on the right-hand of (5.60), resulting to 

r» 777-1 o 



|pliL(t„^i,t„;^) ^ l|ec(0)||^ + 9(51 '*"^" 



77 = 1 



- 2 Xl(^« + Cr7 + /377 + 777 + Sn) Vt^ \\e. 



cllL2(t„-i,t„;J^) 



71 = 1 



+ 2 ^ [{0n + Cn + Pn + In + 5n) V^ + II ^c |1 L2(t„_i,t„;J^)) II PIIl^ (0,t„;J^) ' 

(5.61) 



n=l 
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which is an inequahty of the form 

lal^ <c^ + d-b+{d + b)-a, (5.62) 

where a,b,d ^ IR^+i and c E R arc appropriately chosen. It foUows that 

|a| <max{|c|,|d|} + |d| + |6|, (5.63) 

which, using the notation introduced in the statement of the theorem, imphes 



\p\\h2{o,t,„;y) ^ l|ec(0)|| + W- Y^ K„r„ + 3?7p 



'=c|lL2(0,t„;J^) 



(5.64) 



To close the estimate, the last term on the right-hand side of (5.64) is bounded by 
The first term yields 



L2(0,t„;^) 



E 



|/„_i(«;!^-f/") + Z„e^ 






— Z^ I ln-l£n-l + 'n^n 

^ m ^ 



(5.66) 



n=l 



Similarly, the second term yields 



\\Ua\\ 



L2{0,t^;y) 



< 



_. m 



(5.67) 



Merging these inequalities with (5.641 and using the triangle inequality we obtain 

l|e|lL2(0,t„;J^) < l|e|lL2(0,t„;J^) + IIPllL2(0,t„;^) (5-68) 



we obtain (5.34). 



5.14. Remark (short-time integration). In the spirit of Theorem 3.7 it is pos- 
sible to modify Theorem |5.9| and the appropriate indicators as to accommodate 
a short time-integration version of this result where Li-accumulation in time re- 
places the L2-accumulation for certain estimators. Over shorter time-intervals this 
provides a tighter bound. 

5.15. Theorem (a posteriori energy-error bound for Euler-IPDG). Under the 
same assumptions of Theorem \5.9\ assuming we employ the IPDG method in space 
as described in fpl the error bound {5.34) holds with the estimators r?p,m, 770, m and 
Sr=i ^n'^ri explicitly computable as follows: 

(a) for the elliptic indicators e„ and £+_i, replace S' by Sip, as dehned in {4.9), 
into \5.23) and {5.24}, respectively; 

(b) for the nonconforming part indicators Sn and (5„, respectively, (cf [5.3 and 
Lemma 4.3), we replace 



c/w" - u: 



d+ 



H-i(r2) 



\u. 



d \\a 



and 



u: 



n-l 

d+ 



respectively, by 



CpfCi||v^„|[/"-c/"-i1||, C2||v^ic/"1| 



(5.69) 
(5.70) 



where h„ is the mesh-size function of ^n 9,nd (T„ is related to it via {4.2); 
(c) replace all II~^(i7) norms by Cpp times the L2(^) norm. 
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6. Computer experiments 

In this final section we summarize the resuhs of computer experiments aimed at 
testing the efficiency and reliability of the fully discrete estimators derived in § [5J 
We built our code upon the free finite element software FEniCS [3T] while Matlab® 
was used as an end-tool to visualize the time-behavior of various estimators. 

All the computational examples are in space dimension d ~ 2 and their choice is 
such as to illustrate as many aspect as possible of practical convergence rate (also 
known as experimental order of convergence, in short EOC) and the efFectivity index 
(EI), on uniform space-time meshes, of the proposed a posteriori error indicators 
defined in S 15.41 



6.1. Benchmark solutions. We consider three benchmark problems for which uq 



and / are chosen so that the exact solution u of problem (1.3) coincides with one 
of the following benchmark solutions: 

ui{x,y,t) ~ sin(nt)sm (ttx) sin (ttj/), (6-1) 

U2{x, y, t) = U2{r, 0, t) = sin(^t)(r2 cos2(0) - lf{r^ sin'icj)) - lfr'°g{(l)), (6.2) 



w,3(a;,y, t) = sin(207rt) sin {nx)sin {iry) 



for t G [0, 1] and 



^"'"^M(-l,l)n[0,l)x(-l,0] 



in (6.11 and (6.3) 



in (6.2). 



(6.3) 



(6.4) 



To complete the definition of U2 in (6.2) we consider 
zo := 0.544483736782464 such that sin^{zoUj) = z^sin2(w), with w = 37r/2, (6.5) 
and 



gW ■■= 



sin((zo - l)w) 

1 Zq + 



Zo - i " ' ' Zo + 1 
X (cos((zo - 1)0) - cos((zo + 1)(/))) 
1 1 



sin((zo + l)w) 



-sin((zo-l)(?!)) 



Zo - 1 ' ' " " ' Zo + 1 
X (cos((zo — l)w) — cos((zo + l)w)). 



sin((zo + 1)(/)) 



(6.6) 



It is well-known [22", 'F that the gradient of U2 in (6.2) has a singularity at the 
reentrant corner located at the origin of fl. 

Solution u\ is smooth and varies "slowly" in time. Solution uj, is also smooth by 
it oscillates much faster and is used to emphasize the time-error indicator appearing 
in the parabolic error estimator ?7p.m, defined in (5.35). 



Similar examples have been studied elsewhere, for example in [291 , 130 ) . 
Note that the diffusion tensor, a{x,t\ is a constant function (equal to 1) of 
space-time and that the initial error ||m(0) — t/(0)||^ = in all examples. 

6.2. Computed quantities. In each of the examples, we compute the solution of 



( 2.20 ) using finite element spaces consisting of polynomials of degree -p equal to 1, 2 



and 3 with interior penalty parameter Ca,fj,{3r) in (4.2 ) having values 40, 80 and 160 
respectively which are sufficient to guarantee stability of the numerical scheme. 

We study the asymptotic behavior of the indicators by setting all constants 
appearing in Theorem |5.15| equal to 1 and monitoring the evolution of the values 
and experimental order of convergence of the estimators and the error as well as 
effectivity index over time on a sequence of uniformly refined meshes with a fixed 
time step r and polynomial degree p. For this purpose, we define experimental 
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order of convergence, in symbols EOC, of a given sequence of positive quantities 
a{i) defined on a sequence of meshes of size h{i) by 

vnn( ■\ log(a(» + !)/«(»)) ,^„. 

^^^(«'*) = iog(M* + i)/MO) ^'-'^ 

and the inverse effectivity index, EI, by 

^^^ ML.(0,t^;T) ^ (6.8) 

We use the inverse effectivity index, instead of the (direct) effectivity index, be- 
cause it is easier to visuahze while conveying the same information. It also has the 
advantage of relating directly to the constants appearing in Theorem |5.15[ 

6.3. Conclusions. The numerical experiments clearly indicate that the error esti- 
mators are reliable (as expected from the theory) and efficient. This is clearly seen 
by the match in EOC between the error and the two main estimators rj^^m and ryp.m 
for each m. 

Since we use time-invariant finite element spaces, the mesh-change estimators 
are null and do not infiucnce the estimators. 

The nonconforming indicator {\'Y^^=i^n'''nf^'^sts found to be of higher order 
with respect the elliptic estimator, jye.m ■ This is most likely to be an effect of using 
time-invariant meshes and the nonconforming indicator can be safely ignored as 
long as the mesh does not change. 

Adding mesh change, space-time-dependent diffusion a, and variable time-step 
to our numerical experiments will exhibit more properties of the estimators but 
we eschew deeper numerical experiments in this paper for concision's sake. For the 
same reason, the derivation of adaptive methods based on our indicators is omitted 
here. 



Results for problem ( 6.1 ) with p —1, problem (6.2 ) with p = 2 and problem ( 6.3 ) 



with p = 3 are depicted and commented further in figures [T] [5] and [3] respectively. 
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Figure 1. Example with exact solution ui, given by (6.1), ap- 
proximated with piecewise polynomials of degree p — 1. 
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(a) Mesh-size h{i) = 2-'/2, i = 2, . . . , 8, and timestep t = 0.1 h. On top we plot the EOC 
of the single cumulative indicators r;p,m and rje^m- Both indicators have the same asymptotic 
EOC ft! 1 as has the error. The effectivity index tends towards 1/0.12. 
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(b) Mesh-size h{i) = 2'*/^, i = 2,..., 8, and timestep t = O.l/i^. On top wc plot the EOC 
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Figure 2. Example (6.2) with piecewise polynomials of degree 
p = 2. 
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(a) Mesh-size h{i) = 2~'/-^, i = 2, . . . , 5, and timestep r = 0.1 h'^. On top we plot the EOC 
of the single cumulative indicators r)p,m and rie,m- EOC < 1 for the error is due to lack of 
i/^-rcgularity. Note that the elliptic estimator has asymptotically the same EOC as the error. 
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(b) Mesh-size h(i) = 2'^^'^, i = 2, . . . , 4 and timestep t = O.l/i^. On top we plot the EOC 
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Figure 3. Example (6.3) with discontinuous piecewise polyno- 
mials of degree p = 3. 
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(a) Mesh-size h{i) = 2''/^, i = 2, . . . , 5 and timestep t = 0.1 /i^. On top we plot the EOC 
of the single cumulative indicators r;p,m and rje^m- Both indicators have the same asymptotic 
EOC ft! 1 as has the error. The effectivity index tends towards asymptotic value 200. 
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(b) Mesh-size h{i) = 2-'/^, i = 2, ... ,4 and timestep r = 0.1 h^. On top we plot the EOC 
of the single cumulative indicators ??p,m and rjc^m- Both indicators have the same asymptotic 
EOC Si 1 as has the error. The effectivity index tends towards asymptotic value 200 
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